{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Updates coordinates in Yang precipitation files.  Individual station trajectories are written to a separate files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_yang(fili):\n",
    "    \"\"\"\n",
    "    Reads Excel file containing Yang corrected monthly precipitation\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    fili - file path\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Pandas dataframe containing monthly precipitation\n",
    "    \"\"\"\n",
    "    df = pd.read_excel(fili, sheet_name='monthly-all', header=0, skiprows=[1,2,3], \n",
    "                   na_values='-', usecols=14)\n",
    "    df = df.dropna(how='all')\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_position(fili):\n",
    "    \"\"\"\n",
    "    Reader for position files contained in the NPSNOW dataset\n",
    "\n",
    "    Arguments\n",
    "    ---------\n",
    "    fili - file path\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Pandas dataframe containing drifting station positions\n",
    "    \"\"\"\n",
    "\n",
    "    import calendar\n",
    "    import datetime as dt\n",
    "    \n",
    "    df = pd.read_csv(fili, header=None, delim_whitespace=True,\n",
    "                     names=['year','month','day','hour','lat','lon'])\n",
    "    df['hour'][df['hour'] == 24] = 0 #There is no hour 24\n",
    "    if (df['hour'] > 24).any():\n",
    "        df['hour'][df['hour'] > 24] = 12\n",
    "\n",
    "    # This is a fix to deal with non-valid dates: e.g. 30 February\n",
    "    isday = [row[1]['day'] <= \\\n",
    "             calendar.monthrange( int(row[1]['year']),int(row[1]['month']) )[1] \\\n",
    "             for row in df.iterrows()]\n",
    "    df = df[isday] # only return rows with valid date\n",
    "\n",
    "    df.index = [dt.datetime(int( '19{:2d}'.format(row[1]['year']) ),\n",
    "                            int(row[1]['month']),\n",
    "                            int(row[1]['day']),\n",
    "                            int(row[1]['hour'])) \\\n",
    "                for row in df.iterrows()]\n",
    "    df['lat'] = df['lat'].floordiv(1000).astype(float) + \\\n",
    "                df['lat'].mod(1000).divide(600)\n",
    "    df['lon'] = df['lon'].floordiv(1000).astype(float) + \\\n",
    "                df['lon'].mod(1000).divide(600)\n",
    "    \n",
    "    return df[['lat','lon']]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def llaverage(x):\n",
    "    if (x['lon'].max() > 270.) and (x['lon'].min() < 90.):\n",
    "        xd = x.copy()\n",
    "        xd['lon'].where(xd['lon'] < 180., xd['lon']-360., inplace=True)\n",
    "        result = xd.mean()\n",
    "        if result['lon'] < 0.: result['lon'] = result['lon']+360.\n",
    "        return result\n",
    "    else:\n",
    "        return x.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get Yang precipitation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ediri = r'C:\\Users\\apbarret\\Documents\\data\\Arctic_precip'\n",
    "efile = r'NPP-yang_copy_apb.xls'\n",
    "ppt_df = read_yang(os.path.join(ediri,efile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>NP</th>\n",
       "      <th>YY</th>\n",
       "      <th>MM</th>\n",
       "      <th>ND</th>\n",
       "      <th>Lat</th>\n",
       "      <th>Lon</th>\n",
       "      <th>Tmn</th>\n",
       "      <th>Ug</th>\n",
       "      <th>DP</th>\n",
       "      <th>Dtc</th>\n",
       "      <th>snow%</th>\n",
       "      <th>Pg</th>\n",
       "      <th>windC</th>\n",
       "      <th>traceC</th>\n",
       "      <th>Pc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>31.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>74.51</td>\n",
       "      <td>-150.67</td>\n",
       "      <td>-31.2</td>\n",
       "      <td>4.7</td>\n",
       "      <td>5.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.5</td>\n",
       "      <td>2.5</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>31.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>28.0</td>\n",
       "      <td>74.70</td>\n",
       "      <td>-147.57</td>\n",
       "      <td>-20.3</td>\n",
       "      <td>4.7</td>\n",
       "      <td>19.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>21.3</td>\n",
       "      <td>20.3</td>\n",
       "      <td>0.8</td>\n",
       "      <td>42.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>31.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>75.30</td>\n",
       "      <td>-146.08</td>\n",
       "      <td>-28.8</td>\n",
       "      <td>3.9</td>\n",
       "      <td>6.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>4.3</td>\n",
       "      <td>4.3</td>\n",
       "      <td>2.4</td>\n",
       "      <td>11.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>31.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>75.48</td>\n",
       "      <td>-146.26</td>\n",
       "      <td>-22.2</td>\n",
       "      <td>3.6</td>\n",
       "      <td>9.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>5.1</td>\n",
       "      <td>2.6</td>\n",
       "      <td>2.0</td>\n",
       "      <td>9.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>31.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>75.46</td>\n",
       "      <td>-146.25</td>\n",
       "      <td>-11.2</td>\n",
       "      <td>3.2</td>\n",
       "      <td>18.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>7.4</td>\n",
       "      <td>2.9</td>\n",
       "      <td>1.2</td>\n",
       "      <td>11.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     NP    YY   MM    ND    Lat     Lon   Tmn   Ug    DP   Dtc  snow%    Pg  \\\n",
       "0  31.0  89.0  1.0  31.0  74.51 -150.67 -31.2  4.7   5.0  25.0  100.0   2.0   \n",
       "1  31.0  89.0  2.0  28.0  74.70 -147.57 -20.3  4.7  19.0   8.0  100.0  21.3   \n",
       "2  31.0  89.0  3.0  31.0  75.30 -146.08 -28.8  3.9   6.0  24.0  100.0   4.3   \n",
       "3  31.0  89.0  4.0  30.0  75.48 -146.26 -22.2  3.6   9.0  20.0  100.0   5.1   \n",
       "4  31.0  89.0  5.0  31.0  75.46 -146.25 -11.2  3.2  18.0  12.0  100.0   7.4   \n",
       "\n",
       "   windC  traceC    Pc  \n",
       "0    1.5     2.5   6.0  \n",
       "1   20.3     0.8  42.4  \n",
       "2    4.3     2.4  11.0  \n",
       "3    2.6     2.0   9.7  \n",
       "4    2.9     1.2  11.5  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ppt_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_date(df):\n",
    "    \"\"\"\n",
    "    parses date from Yang DataFrame to be used as index\n",
    "    \"\"\"\n",
    "    date = []\n",
    "    for y, m in zip(df['YY'].astype(int), df['MM'].astype(int)):\n",
    "        if y > 30:\n",
    "            y = y+1900\n",
    "        else:\n",
    "            y = y+2000\n",
    "        date.append(dt.datetime(y,m,1))\n",
    "    return date\n",
    "\n",
    "def extract_station(df, np):\n",
    "    newDf = df[df['NP'] == np].copy()\n",
    "    newDf.index = parse_date(newDf)\n",
    "    newDf.drop(['YY','MM'], axis=1, inplace=True)\n",
    "    return newDf\n",
    "\n",
    "def position_month(np):\n",
    "    pos_diri = r'C:\\Users\\apbarret\\Documents\\data\\SnowOnSeaIce\\NPSNOW'\n",
    "    pos_file = os.path.join(pos_diri,'position','position.{:02d}'.format(np))\n",
    "    position = read_position(pos_file)\n",
    "    mposition = position.resample('MS').agg(llaverage)\n",
    "    return mposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for np in ppt_df.NP.unique().astype(int):\n",
    "    \n",
    "    station = extract_station(ppt_df, np) # Extract Yang data for station NP\n",
    "\n",
    "    if np == 27:\n",
    "        statjn = station # NP27 position file is corrupted so do not update\n",
    "    else:\n",
    "        position = position_month(np) # Get position file and calculate monthly mean position\n",
    "        statjn = station.join(position, rsuffix='_new') # Join DataFrame and Series\n",
    "\n",
    "        # Get coords\n",
    "        x0 = statjn['Lon'].values\n",
    "        y0 = statjn['Lat'].values\n",
    "        x1 = statjn['lon'].values\n",
    "        y1 = statjn['lat'].values\n",
    "\n",
    "        # Replace Lat and Lon\n",
    "        statjn = statjn.drop(['Lat','Lon'], axis=1).rename({'lat': 'Lat', 'lon': 'Lon'}, axis=1) \n",
    "    \n",
    "    # write to file\n",
    "    statjn.to_csv(os.path.join(ediri, 'yang_np_precip_updated_coords_{:02d}.csv'.format(np))) # Write to file\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "    map_proj = ccrs.NorthPolarStereo()\n",
    "\n",
    "    fig = plt.figure(figsize=(10,10))\n",
    "    ax = plt.subplot(projection=map_proj)\n",
    "    ax.set_extent([-180.,180.,65.,90.], ccrs.PlateCarree())\n",
    "    ax.add_feature(cfeature.LAND)\n",
    "    ax.gridlines()\n",
    "\n",
    "    traj0 = map_proj.transform_points(ccrs.PlateCarree(), x0, y0)\n",
    "    traj1 = map_proj.transform_points(ccrs.PlateCarree(), x1, y1)\n",
    "\n",
    "    ax.plot(traj0[:,0], traj0[:,1], '-o', ms=2.5, label='Yang')\n",
    "    ax.plot(traj1[:,0], traj1[:,1], '-o', ms=2.5, label='Recalculated')\n",
    "\n",
    "    ax.legend()\n",
    "    \n",
    "    #fig.savefig(os.path.join(ediri,'yang_np_precip_updated_coords_{:02d}.png'.format(np)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
